ROIs from the Craddock et al. (2012) parcellation atlas

Mean parameter estimates were extracted from parcel 292 and 116.

Extract mean parameter estimates

Run bash script to calculate mean parameter estimates for each subject and condition contrast (condition > rest) within each ROI using AFNI 3dmaskave.

Path to bash script: functional-workshop/code/ROI_analysis/extract_parameterEstimates.sh

#!/bin/bash
. ~/.bashrc

# This script extracts mean parameter estimates and SDs within an ROI or parcel
# from subject FX condition contrasts (condition > rest) for each wave. Output is 
# saved as a text file in the output directory.

# Set paths and variables
# ------------------------------------------------------------------------------------------
# paths

con_dir='/Volumes/psych-cog/dsnlab/MDC/functional-workshop/data/FX_models' #fx contrast directory
atlas_dir='/Volumes/psych-cog/dsnlab/MDC/functional-workshop/data/ROIs' #roi/parcellation atlas directory 
output_dir='/Volumes/psych-cog/dsnlab/MDC/functional-workshop/results/ROI_analysis' #roi/parcellation output directory
rx_model='/Volumes/psych-cog/dsnlab/MDC/functional-workshop/results/AFNI/all+tlrc' #rx model (for atlas alignment only)

# variables
subjects=`cat /Volumes/psych-cog/dsnlab/MDC/functional-workshop/data/subject_list.txt`
parcellation_atlas=(craddock_all.nii.gz) #roi/parcellation atlas file
parcellation_map=(31) #parcellation map number (if applicable)
aligned_parcellation_map=(aligned_craddock_400) #aligned roi/parcellation map name
aligned_parcellation_num=(116 292) #parcellation number(s) to extract from; use $(seq 1 N) where N is the total number of parcels to extract from all
waves=(t1 t2 t3) #waves or task names
fx_cons=(con_0001 con_0002 con_0003 con_0004) #fx con files to extract from

# Align roi/parcellation map to data
# ------------------------------------------------------------------------------------------
echo "aligning parcellation map"
if [ -f $atlas_dir/${aligned_parcellation_map}+tlrc.BRIK ]; then
    echo "aligned parcellation map already exists"
else 
3dAllineate -source $atlas_dir/$parcellation_atlas[$parcellation_map] -master $rx_model -final NN -1Dparam_apply '1D: 12@0'\' -prefix $atlas_dir/$aligned_parcellation_map
fi

# Extract mean parameter estimates and SDs for each subject, wave, contrast, and roi/parcel
# ------------------------------------------------------------------------------------------

for sub in ${subjects[@]}; do 
    for wave in ${waves[@]}; do 
        for con in ${fx_cons[@]}; do 
            for num in ${aligned_parcellation_num[@]}; do 
                echo ${sub} ${wave} ${con} ${num} `3dmaskave -sigma -quiet -mrange $num $num -mask $atlas_dir/${aligned_parcellation_map}+tlrc $con_dir/${sub}_${wave}_${con}.nii` >> $output_dir/parameterEstimates.txt
            done
        done
    done
done

The output will be saved in a text file functional-workshop/results/ROI_analysis/parameterEstimates.txt

library(dplyr)
library(knitr)

read.table('/Volumes/psych-cog/dsnlab/MDC/functional-workshop/results/ROI_analysis/parameterEstimates.txt', sep = "", fill = TRUE, stringsAsFactors=FALSE) %>%
  head(10) %>%
  kable(format = 'pandoc')
V1 V2 V3 V4 V5 V6
s001 t1 con_0001 116 -0.5871270 0.274073
s001 t1 con_0001 292 -0.3261990 0.331749
s001 t1 con_0002 116 -0.0216224 0.342096
s001 t1 con_0002 292 -0.0078532 0.387301
s001 t1 con_0003 116 -0.5910870 0.330697
s001 t1 con_0003 292 -0.2818540 0.396804
s001 t1 con_0004 116 -0.0967547 0.305696
s001 t1 con_0004 292 -0.0509896 0.374470
s001 t2 con_0001 116 -0.6617310 0.311137
s001 t2 con_0001 292 -0.4697610 0.275190

Load packages

library(tidyverse)
library(lme4)
library(lmerTest)
library(wesanderson)

Load data

# load parameter estimate .csv file
data = read.table('/Volumes/psych-cog/dsnlab/MDC/functional-workshop/results/ROI_analysis/parameterEstimates.txt', sep = "", fill = TRUE, stringsAsFactors=FALSE)

# load age covariates and rename variables
age = read.csv('/Volumes/psych-cog/dsnlab/MDC/functional-workshop/data/covariates/age.csv') %>%
  rename("subjectID" = Subj,
         "wave" = wavenum)

Tidy data

Specify your variables names and levels

# tidy raw data
data1 = data %>% 
  # rename variables
  rename('subjectID' = V1,
         'wave' = V2,
         'con' = V3,
         'parcellation' = V4,
         'beta' = V5,
         'sd' = V6) %>%
  # convert con file names to condition names
  mutate(target = ifelse(con %in% c('con_0001', 'con_0002'), 'self', 'other'), 
         domain = ifelse(con %in% c('con_0001', 'con_0003'), 'academic', 'social'), 
  # change data type to factor
         parcellation = as.factor(parcellation),
         target = as.factor(target),
         domain = as.factor(domain)) %>%
  # change to integer
  extract(wave, 'wave', 't([0-3]{1})') %>%
  mutate(wave = as.integer(wave))

Merge data, add age to the data frame and center

merged = left_join(data1, age, by = c('subjectID', 'wave')) %>%
  mutate(age_c = age-mean(age, na.rm=TRUE))

# print data frame header
merged %>%
  head(16) %>%
  kable(format = 'pandoc')
subjectID wave con parcellation beta sd target domain age age_c
s001 1 con_0001 116 -0.5871270 0.274073 self academic 10.41530 -2.0912751
s001 1 con_0001 292 -0.3261990 0.331749 self academic 10.41530 -2.0912751
s001 1 con_0002 116 -0.0216224 0.342096 self social 10.41530 -2.0912751
s001 1 con_0002 292 -0.0078532 0.387301 self social 10.41530 -2.0912751
s001 1 con_0003 116 -0.5910870 0.330697 other academic 10.41530 -2.0912751
s001 1 con_0003 292 -0.2818540 0.396804 other academic 10.41530 -2.0912751
s001 1 con_0004 116 -0.0967547 0.305696 other social 10.41530 -2.0912751
s001 1 con_0004 292 -0.0509896 0.374470 other social 10.41530 -2.0912751
s001 2 con_0001 116 -0.6617310 0.311137 self academic 13.48767 0.9810956
s001 2 con_0001 292 -0.4697610 0.275190 self academic 13.48767 0.9810956
s001 2 con_0002 116 -0.3545030 0.338036 self social 13.48767 0.9810956
s001 2 con_0002 292 -0.1359760 0.424825 self social 13.48767 0.9810956
s001 2 con_0003 116 -0.5791710 0.404919 other academic 13.48767 0.9810956
s001 2 con_0003 292 -0.5360000 0.272914 other academic 13.48767 0.9810956
s001 2 con_0004 116 -0.2961900 0.273353 other social 13.48767 0.9810956
s001 2 con_0004 292 -0.2756520 0.265749 other social 13.48767 0.9810956

Remove missing data to run LME models

data.complete = merged %>%
  na.omit(.)

# print number of rows
cat('raw data: ', nrow(merged))
## raw data:  1944
cat('complete data: ', nrow(data.complete))
## complete data:  1296

Run LME models within parcel 292 and compare

Predict parameter estimates from task conditions (target and domain) and age within parcel 292.

Linear effect of age, random intercepts only

model.1 = lmer(beta ~ target*domain*age_c + (1 | subjectID), data=filter(data.complete, parcellation == 292))
summary(model.1)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: beta ~ target * domain * age_c + (1 | subjectID)
##    Data: filter(data.complete, parcellation == 292)
## 
## REML criterion at convergence: 1106.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0561 -0.4861  0.0324  0.4626  4.3648 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  subjectID (Intercept) 0.1008   0.3175  
##  Residual              0.2582   0.5082  
## Number of obs: 648, groups:  subjectID, 81
## 
## Fixed effects:
##                                Estimate Std. Error        df t value
## (Intercept)                    -0.07709    0.05422 199.40000  -1.422
## targetself                      0.03147    0.05663 554.60000   0.556
## domainsocial                    0.09379    0.05663 554.60000   1.656
## age_c                          -0.03904    0.01658 584.80000  -2.354
## targetself:domainsocial         0.09714    0.08008 554.60000   1.213
## targetself:age_c                0.06418    0.02278 554.60000   2.817
## domainsocial:age_c              0.02990    0.02278 554.60000   1.312
## targetself:domainsocial:age_c  -0.06094    0.03222 554.60000  -1.891
##                               Pr(>|t|)   
## (Intercept)                    0.15660   
## targetself                     0.57861   
## domainsocial                   0.09824 . 
## age_c                          0.01889 * 
## targetself:domainsocial        0.22563   
## targetself:age_c               0.00502 **
## domainsocial:age_c             0.19000   
## targetself:domainsocial:age_c  0.05911 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trgtsl dmnscl age_c  trgts: trgt:_ dmns:_
## targetself  -0.522                                          
## domainsocil -0.522  0.500                                   
## age_c        0.092 -0.052 -0.052                            
## trgtslf:dmn  0.369 -0.707 -0.707  0.037                     
## trgtslf:g_c -0.040  0.076  0.038 -0.687 -0.054              
## domnscl:g_c -0.040  0.038  0.076 -0.687 -0.054  0.500       
## trgtslf:d:_  0.028 -0.054 -0.054  0.486  0.076 -0.707 -0.707

Linear effect of age, random intercepts and age slopes

model.2 = lmer(beta ~ target*domain*age_c + (1 + age_c | subjectID), data=filter(data.complete, parcellation == 292))
summary(model.2)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: beta ~ target * domain * age_c + (1 + age_c | subjectID)
##    Data: filter(data.complete, parcellation == 292)
## 
## REML criterion at convergence: 1101.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8144 -0.4814  0.0320  0.4554  4.3703 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev. Corr 
##  subjectID (Intercept) 0.092839 0.30469       
##            age_c       0.001157 0.03402  -0.60
##  Residual              0.252104 0.50210       
## Number of obs: 648, groups:  subjectID, 81
## 
## Fixed effects:
##                                Estimate Std. Error        df t value
## (Intercept)                    -0.07347    0.05281 204.40000  -1.391
## targetself                      0.03147    0.05595 516.30000   0.562
## domainsocial                    0.09379    0.05595 516.30000   1.676
## age_c                          -0.03807    0.01680 307.30000  -2.266
## targetself:domainsocial         0.09714    0.07913 516.30000   1.228
## targetself:age_c                0.06418    0.02251 516.30000   2.851
## domainsocial:age_c              0.02990    0.02251 516.30000   1.328
## targetself:domainsocial:age_c  -0.06094    0.03184 516.30000  -1.914
##                               Pr(>|t|)   
## (Intercept)                    0.16569   
## targetself                     0.57405   
## domainsocial                   0.09429 . 
## age_c                          0.02412 * 
## targetself:domainsocial        0.22013   
## targetself:age_c               0.00453 **
## domainsocial:age_c             0.18475   
## targetself:domainsocial:age_c  0.05617 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trgtsl dmnscl age_c  trgts: trgt:_ dmns:_
## targetself  -0.530                                          
## domainsocil -0.530  0.500                                   
## age_c       -0.003 -0.051 -0.051                            
## trgtslf:dmn  0.375 -0.707 -0.707  0.036                     
## trgtslf:g_c -0.040  0.076  0.038 -0.670 -0.054              
## domnscl:g_c -0.040  0.038  0.076 -0.670 -0.054  0.500       
## trgtslf:d:_  0.028 -0.054 -0.054  0.474  0.076 -0.707 -0.707

Compare models

model.1: beta ~ target * domain * age_c + (1 | subjectID)

model.2: beta ~ target * domain * age_c + (1 + age_c | subjectID)

anova(model.1, model.2) %>%
  `row.names<-`(c('model.1', 'model.2')) %>%
  kable()
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
model.1 10 1083.649 1128.388 -531.8247 1063.649 NA NA NA
model.2 12 1082.726 1136.413 -529.3631 1058.726 4.923206 2 0.0852981

Adding age as a random effect does not significantly improve the model fit.

Visualize raw data

# set color palette
palette = wes_palette("Zissou", 2, type = "continuous")

Plot fitted curves for parcels 292 and 116

Main effect of target

ggplot(data.complete, aes(x = age, 
                          y = beta, 
                          group = interaction(subjectID, target, domain), 
                          color = target)) +
  geom_point(size = .5, alpha = .1) + 
  geom_line(alpha = .1) + 
  geom_line(aes(group=target), size = 1.5, stat = 'smooth', method = 'lm', formula = y ~ poly(x,2)) + 
  facet_wrap(~parcellation, ncol = 2) +
  geom_hline(yintercept = 0, color = 'gray') +
  scale_color_manual(breaks = c('self', 'other'), values = c(self=palette[2], other=palette[1])) +
  scale_x_continuous(breaks=c(10,13,16)) +
  coord_cartesian(ylim=c(-1,1)) +
  theme_minimal(base_size = 18)

Interaction between target and domain

ggplot(data.complete, aes(x = age, 
                          y = beta, 
                          group = interaction(subjectID, target, domain), 
                          color = target, 
                          linetype = domain)) +
  geom_point(size = .5, alpha = .1) + 
  geom_line(alpha = .1) + 
  geom_line(aes(group=interaction(target,domain)), size = 1.5, stat = 'smooth', method = 'lm', formula = y ~ poly(x,2)) + 
  facet_wrap(~parcellation, ncol = 2) +
  geom_hline(yintercept = 0, color = 'gray')+
  scale_color_manual(breaks = c('self', 'other'), values = c(self=palette[2], other=palette[1]))+
  scale_x_continuous(breaks=c(10,13,16)) +
  coord_cartesian(ylim=c(-1,1)) +
  theme_minimal(base_size = 18)

Plot LOESS curves for parcels 292 and 116

Main effect of target

ggplot(data.complete, aes(x = age, 
                          y = beta, 
                          group = interaction(subjectID, target, domain), 
                          color = target, 
                          linetype = domain)) +
  geom_point(size = .5, alpha = .1) + 
  geom_line(alpha = .1) + 
  geom_line(aes(group=target), size = 1.5, stat = 'smooth', method = 'loess') + 
  facet_wrap(~parcellation, ncol = 2) +
  geom_hline(yintercept = 0, color = 'gray')+
  scale_color_manual(breaks = c('self', 'other'), values = c(self=palette[2], other=palette[1]))+
  scale_x_continuous(breaks=c(10,13,16)) +
  coord_cartesian(ylim=c(-1,1)) +
  theme_minimal(base_size = 18)

Interaction between target and domain

ggplot(data.complete, aes(x = age, 
                          y = beta, 
                          group = interaction(subjectID, target, domain), 
                          color = target, 
                          linetype = domain)) +
  geom_point(size = .5, alpha = .1) + 
  geom_line(alpha = .1) + 
  geom_line(aes(group=interaction(target,domain)), size = 1.5, stat = 'smooth', method = 'loess') + 
  facet_wrap(~parcellation, ncol = 2) +
  geom_hline(yintercept = 0, color = 'gray')+
  scale_color_manual(breaks = c('self', 'other'), values = c(self=palette[2], other=palette[1]))+
  scale_x_continuous(breaks=c(10,13,16)) +
  coord_cartesian(ylim=c(-1,1)) +
  theme_minimal(base_size = 18)

Visualize predicted values from model

Plot fitted curves for parcels 292 and 116

# extract random effects formula from model.2 and reconstruct it to use with the `predict` function
REFormulaString = as.character(findbars(model.2@call$formula)[[1]])
REFormula = as.formula(paste0('~(', REFormulaString[[2]], REFormulaString[[1]], REFormulaString[[3]], ')'))

# get expected values for each observation based on model.2
data.complete$expected <- predict(model.2, newdata = data.complete, re.form=REFormula)
data.complete$expected_mean <- predict(model.2, newdata = data.complete, re.form=NA)

Main effect of target

ggplot(data.complete, aes(x = age, 
                          y = expected_mean, 
                          color = target)) +
  geom_line(size = 1.5, stat = 'smooth', method = 'lm', formula = y ~ poly(x,2)) + 
  facet_wrap(~parcellation, ncol = 2) +
  geom_hline(yintercept = 0, color = 'gray')+
  scale_color_manual(breaks = c('self', 'other'), values = c(self=palette[2], other=palette[1]))+
  scale_x_continuous(breaks=c(10,13,16)) +
  coord_cartesian(ylim=c(-1,1)) +
  theme_minimal(base_size = 18)

Interaction between target and domain

ggplot(data.complete, aes(x = age, 
                          y = expected_mean, 
                          group = interaction(target, domain), 
                          color = target, 
                          linetype = domain)) +
  geom_line(size = 1.5, stat = 'smooth', method = 'lm', formula = y ~ poly(x,2)) + 
  facet_wrap(~parcellation, ncol = 2) +
  geom_hline(yintercept = 0, color = 'gray')+
  scale_color_manual(breaks = c('self', 'other'), values = c(self=palette[2], other=palette[1]))+
  scale_x_continuous(breaks=c(10,13,16)) +
  coord_cartesian(ylim=c(-1,1)) +
  theme_minimal(base_size = 18)